subroutine qvflux_interp99(nsub,eallz,eally,eallx,eallzs,eallys,eallxs,eallz_new,eally_new,eallx_new,&
    eallz_news,eally_news,eallx_news,u,v,w,qv,tq,rho,level,dx,dy,wtrhold,qtrhold,&
    qvinflux_x_iz,qvinflux_y_iz,qvinflux_z_iz,qvinflux_iz,qvoutflux_x_iz,qvoutflux_y_iz,qvoutflux_z_iz,qvoutflux_iz,&
    qvflux_e_iz,qvflux_d_iz,v_iz)
    implicit none
    integer, intent(in) :: nsub
    integer, intent(in) :: eallz,eally,eallx,eallzs,eallys,eallxs
    integer, intent(in) :: eallz_new,eally_new,eallx_new,eallz_news,eally_news,eallx_news
    
    real, intent(in) :: dx
    real, intent(in) :: dy
    real, intent(in) :: wtrhold
    real, intent(in) :: qtrhold
    real, intent(in) :: u  (eallx, eallys, eallzs)
    real, intent(in) :: v  (eallxs, eally, eallzs)
    real, intent(in) :: w  (eallxs, eallys, eallz)
    real, intent(in) :: qv   (eallx, eally, eallz)
    real, intent(in) :: tq(eallxs, eallys, eallzs)
    real, intent(in) :: rho               (eallzs)
    real, intent(in) :: level              (eallz)
    
    real, intent(out) :: qvinflux_x_iz(eallzs)
    real, intent(out) :: qvinflux_y_iz(eallzs)
    real, intent(out) :: qvinflux_z_iz(eallzs)
    real, intent(out) :: qvinflux_iz(eallzs)
    real, intent(out) :: qvoutflux_x_iz(eallzs)
    real, intent(out) :: qvoutflux_y_iz(eallzs)
    real, intent(out) :: qvoutflux_z_iz(eallzs)
    real, intent(out) :: qvoutflux_iz(eallzs)
    real, intent(out) :: qvflux_e_iz(eallzs)
    real, intent(out) :: qvflux_d_iz(eallzs)
    real, intent(out) :: v_iz(eallzs)
    
    ! Local variables@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
    real :: qvsx(eallx, eallys, eallzs)
    real :: qvsy(eallxs, eally, eallzs)
    real :: qvsz(eallxs, eallys, eallz)
    real :: qvvx(eallx, eallys, eallzs)
    real :: qvvy(eallxs, eally, eallzs)
    real :: qvvz(eallxs, eallys, eallz)
    real :: u_new(eallx_new, eally_news, eallz_news)
    real :: v_new(eallx_news, eally_new, eallz_news)
    real :: w_new(eallx_news, eally_news, eallz_new)
    real :: qvvx_new(eallx_new, eally_news, eallz_news)
    real :: qvvy_new(eallx_news, eally_new, eallz_news)
    real :: qvvz_new(eallx_news, eally_news, eallz_new)
    
    real :: ws_new(eallx_news, eally_news, eallz_new)
    real :: tq_new(eallx_news, eally_news, eallz_news)
    real :: Wx(eallx_news, eally_news, eallz_news)
    real :: Wy(eallx_news, eally_news, eallz_news)
    real :: Wz(eallx_news, eally_news, eallz_news)
    logical :: cloud_mask(eallx_news, eally_news, eallz_news)
    real :: qvinflux_x_iz_new(eallz_news)
    real :: qvinflux_y_iz_new(eallz_news)
    real :: qvinflux_z_iz_new(eallz_news)
    real :: qvinflux_iz_new(eallz_news)
    real :: qvoutflux_x_iz_new(eallz_news)
    real :: qvoutflux_y_iz_new(eallz_news)
    real :: qvoutflux_z_iz_new(eallz_news)
    real :: qvoutflux_iz_new(eallz_news)
    real :: qvflux_e_iz_new(eallz_news)
    real :: qvflux_d_iz_new(eallz_news)
    real :: v_iz_new(eallz_news)
    real :: flux
    real :: influx_east
    real :: outflux_east
    real :: influx_west
    real :: outflux_west
    real :: influx_north
    real :: outflux_north
    real :: influx_south
    real :: outflux_south
    real :: influx_top
    real :: outflux_top
    real :: influx_bottom
    real :: outflux_bottom
    real :: qvflux_diff
    
    real :: height(eallzs)
    real :: height_new(eallz_news)
    real :: rho_new(eallz_news)
    real :: rho_newz(eallz_new)
    
    integer :: iz,iy,ix,i,irs,ire
    integer :: numiz
    
    !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
    !
    ! Beginning of the qvflux_interpolate code
    !
    !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

    qvsx = qv(:, 1:eallys, 1:eallzs)
    qvsy = qv(1:eallxs, :, 1:eallzs)
    qvsz = qv(1:eallxs, 1:eallys, :)
    qvvx = 0.0
    qvvy = 0.0
    qvvz = 0.0
    
    ! Interpolate Qv to vector point at parent grid
    !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
    call avgx(qvsx, 1, eallx,eallys,eallzs, 2,eallx,1,eallys,1,eallzs, qvvx)
    call avgy(qvsy, 1, eallxs,eally,eallzs, 1,eallxs,2,eally,1,eallzs, qvvy)
    call avgz(qvsz, 1, eallxs,eallys,eallz, 1,eallxs,1,eallys,2,eallz, qvvz)
    qvvx(1,:,:) = qvvx(2,:,:)
    qvvy(:,1,:) = qvvy(:,2,:)
    qvvz(:,:,1) = qvvz(:,:,2)
    
    ! interpolate wind and qv to nest grid
    !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
    call interpolate_uqv(u, qvvx, u_new, qvvx_new, eallxs, eallys, eallzs, nsub)
    call interpolate_vqv(v, qvvy, v_new, qvvy_new, eallxs, eallys, eallzs, nsub)
    call interpolate_wqv(w, qvvz, w_new, qvvz_new, eallxs, eallys, eallzs, nsub)
    ! For test###################################################
    !print*, 'maxval of w_new'
    !write(*,'(3F10.3)') maxval(w_new)
    !print*, 'maxval of u_new'
    !write(*,'(3F10.3)') maxval(u_new)
    !print*, 'maxval of v_new'
    !write(*,'(3F10.3)') maxval(v_new)
    
    ! interpolate w and tq to scalar point at nest grid
    !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
    ws_new = 0.0
    call avgz(w_new, 0, eallx_news,eally_news,eallz_new, 1,eallx_news,1,eally_news,1,eallz_new-1,ws_new)
    ws_new(:,:,eallz_new) = ws_new(:,:,eallz_new-1)
    call interpolate_qvs(tq, tq_new, eallxs, eallys, eallzs, nsub)
    
    ! Deal with rho and Wx,Wy,Wz
    !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
    do iz = 1,eallzs
        height(iz) = level(iz+1)-level(iz)
        irs = (iz-1)*nsub+1
        ire = iz*nsub
        rho_new(irs:ire) = rho(iz)
        Wx(:,:,irs:ire) = (height(iz) / real(nsub)) * (dx / real(nsub))
        Wy(:,:,irs:ire) = (height(iz) / real(nsub)) * (dy / real(nsub))
        Wz(:,:,irs:ire) = (dx / real(nsub)) * (dy / real(nsub))
        height_new(irs:ire) = height(iz) / real(nsub)
    end do
    call avgz_1d(rho_new, 1, eallz_news, 2,eallz_news, rho_newz)
    rho_newz(1) = rho_newz(2)
    rho_newz(eallz_new) = rho_newz(eallz_news) ! Only need to interpolate rho_new at 'Z' direction.
                                               ! It's homogeneous horizontally.
    
    
    !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
    !                                                                                    @
    ! Calculate qv flux at nest grid                                                     @
    !                                                                                    @
    !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
    qvinflux_x_iz_new = 0.0
    qvinflux_y_iz_new = 0.0
    qvinflux_z_iz_new = 0.0
    qvinflux_iz_new = 0.0
    qvinflux_x_iz = 0.0
    qvinflux_y_iz = 0.0
    qvinflux_z_iz = 0.0
    qvinflux_iz = 0.0
    qvoutflux_x_iz_new = 0.0
    qvoutflux_y_iz_new = 0.0
    qvoutflux_z_iz_new = 0.0
    qvoutflux_iz_new = 0.0
    qvoutflux_x_iz = 0.0
    qvoutflux_y_iz = 0.0
    qvoutflux_z_iz = 0.0
    qvoutflux_iz = 0.0
    qvflux_e_iz_new = 0.0
    qvflux_d_iz_new = 0.0
    qvflux_e_iz = 0.0
    qvflux_d_iz = 0.0
    v_iz_new = 0.0
    v_iz = 0.0
    cloud_mask = .false.
    
    !#############################################
    ! Calculate cloud mask
    !#############################################
    do iz = 1,eallz_news
        do iy = 1, eally_news
            do ix = 1, eallx_news
                if ( ws_new(ix,iy,iz) >= wtrhold .and. tq_new(ix,iy,iz) >= qtrhold ) then
                    cloud_mask(ix,iy,iz) = .true.
                    v_iz_new(iz) = v_iz_new(iz) + height_new(iz) * (dx / real(nsub)) * (dy / real(nsub))
                end if
            end do
        end do
    end do
    
    !#####################################################################################
    !                                                                                    #
    ! Start Calculate Qv Flux through Cloud Boundary                                     #
    !                                                                                    #
    !#####################################################################################
    do iz = 1,eallz_news
        do iy = 1, eally_news
            do ix = 1, eallx_news
                    
                !#####################################################
                ! Check if neighboeing girds belong to cloud
                !#####################################################
                if (.not. cloud_mask(ix,iy,iz)) cycle  ! Skep non-cloud cell
                    
                    outflux_east = 0.0
                    influx_east = 0.0
                    outflux_west = 0.0
                    influx_west = 0.0
                    outflux_north = 0.0
                    influx_north = 0.0
                    outflux_south = 0.0
                    influx_south = 0.0
                    outflux_top = 0.0
                    influx_top = 0.0
                    outflux_bottom = 0.0
                    influx_bottom = 0.0
                    qvflux_diff = 0.0
                    flux = 0.0

                    ! --- Check East Boundary ---
                    if (ix < eallx_news) then
                        if (.not. cloud_mask(ix+1,iy,iz)) then
                        flux = u_new(ix+1,iy,iz) * rho_new(iz) * qvvx_new(ix+1,iy,iz) * Wx(ix,iy,iz)  ! Flux = V * rho * Qv * A
                        !qvflux_iz_new(iz) = qvflux_iz_new(iz) + flux
                        !qvflux_x_iz_new(iz) = qvflux_x_iz_new(iz) + flux
                        
                        ! For test###################################################
                        !if (ix <= 300 .and. iy <= 300) then
                        !    write(*,'(3F10.3)') flux,Wx(ix,iy,iz),rho_new(iz)
                        !endif
                          
                            
                        if (flux > 0) then
                            qvoutflux_iz_new(iz) = qvoutflux_iz_new(iz) + flux
                            qvoutflux_x_iz_new(iz) = qvoutflux_x_iz_new(iz) + flux
                            outflux_east = flux
                        endif
                        if (flux < 0) then
                            qvinflux_iz_new(iz) = qvinflux_iz_new(iz) - flux
                            qvinflux_x_iz_new(iz) = qvinflux_x_iz_new(iz) - flux
                            influx_east = (-1.0) * flux
                        endif
                            
                        endif
                    endif

                    ! --- Check West Boundary ---
                    if (ix > 1) then
                        if (.not. cloud_mask(ix-1,iy,iz)) then
                        flux = u_new(ix,iy,iz) * rho_new(iz) * qvvx_new(ix,iy,iz) * Wx(ix,iy,iz)
                            
                        if (flux < 0) then
                            qvoutflux_iz_new(iz) = qvoutflux_iz_new(iz) - flux
                            qvoutflux_x_iz_new(iz) = qvoutflux_x_iz_new(iz) - flux
                            outflux_west = (-1.0) * flux
                        end if
                        if (flux > 0) then
                            qvinflux_iz_new(iz) = qvinflux_iz_new(iz) + flux
                            qvinflux_x_iz_new(iz) = qvinflux_x_iz_new(iz) + flux
                            influx_west = flux
                        end if
                            
                        endif
                    endif

                    ! --- Check North Boundary ---
                    if (iy < eally_news) then
                        if (.not. cloud_mask(ix,iy+1,iz)) then
                        flux = v_new(ix,iy+1,iz) * rho_new(iz) * qvvy_new(ix,iy+1,iz) * Wy(ix,iy,iz)
                            
                        if (flux > 0) then
                            qvoutflux_iz_new(iz) = qvoutflux_iz_new(iz) + flux
                            qvoutflux_y_iz_new(iz) = qvoutflux_y_iz_new(iz) + flux
                            outflux_north = flux
                        endif
                        if (flux < 0) then
                            qvinflux_iz_new(iz) = qvinflux_iz_new(iz) - flux
                            qvinflux_y_iz_new(iz) = qvinflux_y_iz_new(iz) - flux
                            influx_north = (-1.0) * flux
                        endif
                        
                        endif
                    endif

                    ! --- Check South Boundary ---
                    if (iy > 1) then
                        if (.not. cloud_mask(ix,iy-1,iz)) then
                        flux = v_new(ix,iy,iz) * rho_new(iz) * qvvy_new(ix,iy,iz) * Wy(ix,iy,iz)
                            
                        if (flux < 0) then
                            qvoutflux_iz_new(iz) = qvoutflux_iz_new(iz) - flux
                            qvoutflux_y_iz_new(iz) = qvoutflux_y_iz_new(iz) - flux
                            outflux_south = (-1.0) * flux
                        endif
                        if (flux > 0) then
                            qvinflux_iz_new(iz) = qvinflux_iz_new(iz) + flux
                            qvinflux_y_iz_new(iz) = qvinflux_y_iz_new(iz) + flux
                            influx_south = flux
                        endif
                            
                        endif
                    endif

                    ! --- Check Top Boundary ---
                    if (iz < eallz_news) then
                        if (.not. cloud_mask(ix,iy,iz+1)) then
                        flux = w_new(ix,iy,iz+1) * rho_newz(iz+1) * qvvz_new(ix,iy,iz+1) * Wz(ix,iy,iz)
                            
                        if (flux > 0) then
                            qvoutflux_iz_new(iz) = qvoutflux_iz_new(iz) + flux
                            qvoutflux_z_iz_new(iz) = qvoutflux_z_iz_new(iz) + flux
                            outflux_top = flux
                        endif
                        if (flux < 0) then
                            qvinflux_iz_new(iz) = qvinflux_iz_new(iz) - flux
                            qvinflux_z_iz_new(iz) = qvinflux_z_iz_new(iz) - flux
                            influx_top = (-1.0) * flux
                        endif
                            
                        endif
                    endif

                    ! --- Check Bottom Boundary ---
                    if (iz > 1) then
                        if (.not. cloud_mask(ix,iy,iz-1)) then
                        flux = w_new(ix,iy,iz) * rho_newz(iz) * qvvz_new(ix,iy,iz) * Wz(ix,iy,iz)
                            
                        if (flux < 0) then
                            qvoutflux_iz_new(iz) = qvoutflux_iz_new(iz) - flux
                            qvoutflux_z_iz_new(iz) = qvoutflux_z_iz_new(iz) - flux
                            outflux_bottom = (-1.0) * flux
                        endif
                        if (flux > 0) then
                            qvinflux_iz_new(iz) = qvinflux_iz_new(iz) + flux
                            qvinflux_z_iz_new(iz) = qvinflux_z_iz_new(iz) + flux
                            influx_bottom = flux
                        endif
                            
                        endif
                    endif
                        
                    qvflux_diff = (influx_east+influx_west+influx_north+influx_south+influx_top+influx_bottom)-&
                                    (outflux_east+outflux_west+outflux_north+outflux_south+outflux_top+outflux_bottom)
                    if (qvflux_diff >= 0) then
                        qvflux_e_iz_new(iz) = qvflux_e_iz_new(iz) + qvflux_diff
                    else
                        qvflux_d_iz_new(iz) = qvflux_d_iz_new(iz) - qvflux_diff
                    endif

            end do
        end do
    end do
    
    !#############################################
    ! Calculate qvflux at parent grid
    !#############################################
    do iz = 1,eallzs
        irs = (iz-1)*nsub+1
        ire = iz*nsub
        qvinflux_x_iz(iz) = sum(qvinflux_x_iz_new(irs:ire))
        qvinflux_y_iz(iz) = sum(qvinflux_y_iz_new(irs:ire))
        qvinflux_z_iz(iz) = sum(qvinflux_z_iz_new(irs:ire))
        qvinflux_iz(iz) = sum(qvinflux_iz_new(irs:ire))
        
        qvoutflux_x_iz(iz) = sum(qvoutflux_x_iz_new(irs:ire))
        qvoutflux_y_iz(iz) = sum(qvoutflux_y_iz_new(irs:ire))
        qvoutflux_z_iz(iz) = sum(qvoutflux_z_iz_new(irs:ire))
        qvoutflux_iz(iz) = sum(qvoutflux_iz_new(irs:ire))
        
        qvflux_e_iz(iz) = sum(qvflux_e_iz_new(irs:ire))
        qvflux_d_iz(iz) = sum(qvflux_d_iz_new(irs:ire))
        
        v_iz(iz) = sum(v_iz_new(irs:ire))
    end do

    return
end subroutine qvflux_interp99 !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
    

    
    
!############################################################
!                                                           #
! Interpolate U                                             #
!                                                           #
!############################################################       
subroutine interpolate_u(u, u_new, nx, ny, nz, nsub)
  implicit none
  integer, intent(in) :: nx, ny, nz
  integer, intent(in) :: nsub
  real, intent(in) :: u(nx+1, ny, nz)
  real, intent(out) :: u_new(nx*nsub + 1, ny*nsub, nz*nsub)
  real, allocatable :: temp_x(:,:,:), temp_xy(:,:,:)
  integer :: ix, iy, iz, ix_new, iy_new, iz_new, i_old, j_old, k_old
  real :: x_new, y_new, z_new, delta

  allocate(temp_x(nx*nsub + 1, ny, nz))
  allocate(temp_xy(nx*nsub + 1, ny*nsub, nz))

  ! X����߽����ֵ
  do iz = 1, nz
    do iy = 1, ny
      do ix_new = 1, nx*nsub + 1
        x_new = real(ix_new - 1) / real(nsub)
        i_old = floor(x_new) + 1
        if (i_old < 1) i_old = 1
        if (i_old > nx) i_old = nx
        delta = x_new - (i_old - 1)
        if (i_old <= nx) then
          temp_x(ix_new, iy, iz) = u(i_old, iy, iz)*(1.0 - delta) + u(i_old+1, iy, iz)*delta
        else
          temp_x(ix_new, iy, iz) = u(i_old, iy, iz)
        endif
      enddo
    enddo
  enddo

  ! Y�������ĵ��ֵ
  do iz = 1, nz
    do ix_new = 1, nx*nsub + 1
      do iy_new = 1, ny*nsub
        y_new = real(iy_new - 0.5) / real(nsub)
        j_old = floor(y_new + 0.5)
        if (j_old < 1) j_old = 1
        if (j_old > ny) j_old = ny
        delta = y_new - (j_old - 0.5)
        if (j_old < ny) then
          temp_xy(ix_new, iy_new, iz) = temp_x(ix_new, j_old, iz)*(1.0 - delta) + temp_x(ix_new, j_old+1, iz)*delta
        else
          temp_xy(ix_new, iy_new, iz) = temp_x(ix_new, j_old, iz)
        endif
      enddo
    enddo
  enddo

  ! Z�������ĵ��ֵ
  do iy_new = 1, ny*nsub
    do ix_new = 1, nx*nsub + 1
      do iz_new = 1, nz*nsub
        z_new = real(iz_new - 0.5) / real(nsub)
        k_old = floor(z_new + 0.5)
        if (k_old < 1) k_old = 1
        if (k_old > nz) k_old = nz
        delta = z_new - (k_old - 0.5)
        if (k_old < nz) then
          u_new(ix_new, iy_new, iz_new) = temp_xy(ix_new, iy_new, k_old)*(1.0 - delta) + temp_xy(ix_new, iy_new, k_old+1)*delta
        else
          u_new(ix_new, iy_new, iz_new) = temp_xy(ix_new, iy_new, k_old)
        endif
      enddo
    enddo
  enddo

  deallocate(temp_x, temp_xy)
  return
end subroutine interpolate_u !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@



!############################################################
!                                                           #
! Interpolate U and Qv                                      #
!                                                           #
!############################################################       
subroutine interpolate_uqv(u, qv, u_new, qv_new, nx, ny, nz, nsub)
  implicit none
  integer, intent(in) :: nx, ny, nz
  integer, intent(in) :: nsub
  real, intent(in) :: u(nx+1, ny, nz)
  real, intent(in) :: qv(nx+1, ny, nz)
  real, intent(out) :: u_new(nx*nsub + 1, ny*nsub, nz*nsub)
  real, intent(out) :: qv_new(nx*nsub + 1, ny*nsub, nz*nsub)
  real, allocatable :: temp_x(:,:,:), temp_xy(:,:,:), temp_qvx(:,:,:), temp_qvxy(:,:,:)
  integer :: ix, iy, iz, ix_new, iy_new, iz_new, i_old, j_old, k_old
  real :: x_new, y_new, z_new, delta

  allocate(temp_x(nx*nsub + 1, ny, nz))
  allocate(temp_xy(nx*nsub + 1, ny*nsub, nz))
  allocate(temp_qvx(nx*nsub + 1, ny, nz))
  allocate(temp_qvxy(nx*nsub + 1, ny*nsub, nz))
  
  !write(*, '(I5)') nx*nsub,ny*nsub,nz*nsub

  ! X����߽����ֵ
  do iz = 1, nz
    do iy = 1, ny
      do ix_new = 1, nx*nsub + 1
        x_new = real(ix_new - 1) / real(nsub)
        i_old = floor(x_new) + 1
        if (i_old < 1) i_old = 1
        if (i_old > nx) i_old = nx
        delta = x_new - (i_old - 1)
        if (i_old <= nx) then
          temp_x(ix_new, iy, iz) = u(i_old, iy, iz)*(1.0 - delta) + u(i_old+1, iy, iz)*delta
          temp_qvx(ix_new, iy, iz) = qv(i_old, iy, iz)*(1.0 - delta) + qv(i_old+1, iy, iz)*delta
        else
          temp_x(ix_new, iy, iz) = u(i_old, iy, iz)
          temp_qvx(ix_new, iy, iz) = qv(i_old, iy, iz)
        endif
      enddo
    enddo
  enddo

  ! Y�������ĵ��ֵ
  do iz = 1, nz
    do ix_new = 1, nx*nsub + 1
      do iy_new = 1, ny*nsub
        y_new = real(iy_new - 0.5) / real(nsub)
        j_old = floor(y_new + 0.5)
        if (j_old < 1) j_old = 1
        if (j_old > ny) j_old = ny
        delta = y_new - (j_old - 0.5)
        if (j_old < ny) then
          temp_xy(ix_new, iy_new, iz) = temp_x(ix_new, j_old, iz)*(1.0 - delta) + temp_x(ix_new, j_old+1, iz)*delta
          temp_qvxy(ix_new, iy_new, iz) = temp_qvx(ix_new, j_old, iz)*(1.0 - delta) + temp_qvx(ix_new, j_old+1, iz)*delta
        else
          temp_xy(ix_new, iy_new, iz) = temp_x(ix_new, j_old, iz)
          temp_qvxy(ix_new, iy_new, iz) = temp_qvx(ix_new, j_old, iz)
        endif
      enddo
    enddo
  enddo

  ! Z�������ĵ��ֵ
  do iy_new = 1, ny*nsub
    do ix_new = 1, nx*nsub + 1
      do iz_new = 1, nz*nsub
        z_new = real(iz_new - 0.5) / real(nsub)
        k_old = floor(z_new + 0.5)
        if (k_old < 1) k_old = 1
        if (k_old > nz) k_old = nz
        delta = z_new - (k_old - 0.5)
        if (k_old < nz) then
          u_new(ix_new, iy_new, iz_new) = temp_xy(ix_new, iy_new, k_old)*(1.0 - delta) + temp_xy(ix_new, iy_new, k_old+1)*delta
          qv_new(ix_new, iy_new, iz_new) = temp_qvxy(ix_new, iy_new, k_old)*(1.0 - delta) + temp_qvxy(ix_new, iy_new, k_old+1)*delta
        else
          u_new(ix_new, iy_new, iz_new) = temp_xy(ix_new, iy_new, k_old)
          qv_new(ix_new, iy_new, iz_new) = temp_qvxy(ix_new, iy_new, k_old)
        endif
      enddo
    enddo
  enddo

  deallocate(temp_x, temp_xy, temp_qvx, temp_qvxy)
  return
end subroutine interpolate_uqv !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@



!############################################################
!                                                           #
! Interpolate V                                             #
!                                                           #
!############################################################     
subroutine interpolate_v(v, v_new, nx, ny, nz, nsub)
  implicit none
  integer, intent(in) :: nx, ny, nz
  integer, intent(in) :: nsub
  real, intent(in) :: v(nx, ny+1, nz)
  real, intent(out) :: v_new(nx*nsub, ny*nsub + 1, nz*nsub)
  real, allocatable :: temp_y(:,:,:), temp_xy(:,:,:)
  integer :: ix, iy, iz, ix_new, iy_new, iz_new, i_old, j_old, k_old
  real :: x_new, y_new, z_new, delta

  allocate(temp_y(nx, ny*nsub + 1, nz))
  allocate(temp_xy(nx*nsub, ny*nsub + 1, nz))

  ! Y����߽����ֵ
  do iz = 1, nz
    do ix = 1, nx
      do iy_new = 1, ny*nsub + 1
        y_new = real(iy_new - 1) / real(nsub)
        j_old = floor(y_new) + 1
        if (j_old < 1) j_old = 1
        if (j_old > ny) j_old = ny
        delta = y_new - (j_old - 1)
        if (j_old <= ny) then
          temp_y(ix, iy_new, iz) = v(ix, j_old, iz)*(1.0 - delta) + v(ix, j_old+1, iz)*delta
        else
          temp_y(ix, iy_new, iz) = v(ix, j_old, iz)
        endif
      enddo
    enddo
  enddo

  ! X�������ĵ��ֵ
  do iz = 1, nz
    do iy_new = 1, ny*nsub + 1
      do ix_new = 1, nx*nsub
        x_new = real(ix_new - 0.5) / real(nsub)
        i_old = floor(x_new + 0.5)
        if (i_old < 1) i_old = 1
        if (i_old > nx) i_old = nx
        delta = x_new - (i_old - 0.5)
        if (i_old < nx) then
          temp_xy(ix_new, iy_new, iz) = temp_y(i_old, iy_new, iz)*(1.0 - delta) + temp_y(i_old+1, iy_new, iz)*delta
        else
          temp_xy(ix_new, iy_new, iz) = temp_y(i_old, iy_new, iz)
        endif
      enddo
    enddo
  enddo

  ! Z�������ĵ��ֵ
  do iy_new = 1, ny*nsub + 1
    do ix_new = 1, nx*nsub
      do iz_new = 1, nz*nsub
        z_new = real(iz_new - 0.5) / real(nsub)
        k_old = floor(z_new + 0.5)
        if (k_old < 1) k_old = 1
        if (k_old > nz) k_old = nz
        delta = z_new - (k_old - 0.5)
        if (k_old < nz) then
          v_new(ix_new, iy_new, iz_new) = temp_xy(ix_new, iy_new, k_old)*(1.0 - delta) + temp_xy(ix_new, iy_new, k_old+1)*delta
        else
          v_new(ix_new, iy_new, iz_new) = temp_xy(ix_new, iy_new, k_old)
        endif
      enddo
    enddo
  enddo

  deallocate(temp_y, temp_xy)
  return
end subroutine interpolate_v !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@



!############################################################
!                                                           #
! Interpolate V and Qv                                      #
!                                                           #
!############################################################     
subroutine interpolate_vqv(v, qv, v_new, qv_new, nx, ny, nz, nsub)
  implicit none
  integer, intent(in) :: nx, ny, nz
  integer, intent(in) :: nsub
  real, intent(in) :: v(nx, ny+1, nz)
  real, intent(in) :: qv(nx, ny+1, nz)
  real, intent(out) :: v_new(nx*nsub, ny*nsub + 1, nz*nsub)
  real, intent(out) :: qv_new(nx*nsub, ny*nsub + 1, nz*nsub)
  real, allocatable :: temp_y(:,:,:), temp_xy(:,:,:), temp_qvy(:,:,:), temp_qvxy(:,:,:)
  integer :: ix, iy, iz, ix_new, iy_new, iz_new, i_old, j_old, k_old
  real :: x_new, y_new, z_new, delta

  allocate(temp_y(nx, ny*nsub + 1, nz))
  allocate(temp_xy(nx*nsub, ny*nsub + 1, nz))
  allocate(temp_qvy(nx, ny*nsub + 1, nz))
  allocate(temp_qvxy(nx*nsub, ny*nsub + 1, nz))

  ! Y����߽����ֵ
  do iz = 1, nz
    do ix = 1, nx
      do iy_new = 1, ny*nsub + 1
        y_new = real(iy_new - 1) / real(nsub)
        j_old = floor(y_new) + 1
        if (j_old < 1) j_old = 1
        if (j_old > ny) j_old = ny
        delta = y_new - (j_old - 1)
        if (j_old <= ny) then
          temp_y(ix, iy_new, iz) = v(ix, j_old, iz)*(1.0 - delta) + v(ix, j_old+1, iz)*delta
          temp_qvy(ix, iy_new, iz) = qv(ix, j_old, iz)*(1.0 - delta) + qv(ix, j_old+1, iz)*delta
        else
          temp_y(ix, iy_new, iz) = v(ix, j_old, iz)
          temp_qvy(ix, iy_new, iz) = qv(ix, j_old, iz)
        endif
      enddo
    enddo
  enddo

  ! X�������ĵ��ֵ
  do iz = 1, nz
    do iy_new = 1, ny*nsub + 1
      do ix_new = 1, nx*nsub
        x_new = real(ix_new - 0.5) / real(nsub)
        i_old = floor(x_new + 0.5)
        if (i_old < 1) i_old = 1
        if (i_old > nx) i_old = nx
        delta = x_new - (i_old - 0.5)
        if (i_old < nx) then
          temp_xy(ix_new, iy_new, iz) = temp_y(i_old, iy_new, iz)*(1.0 - delta) + temp_y(i_old+1, iy_new, iz)*delta
          temp_qvxy(ix_new, iy_new, iz) = temp_qvy(i_old, iy_new, iz)*(1.0 - delta) + temp_qvy(i_old+1, iy_new, iz)*delta
        else
          temp_xy(ix_new, iy_new, iz) = temp_y(i_old, iy_new, iz)
          temp_qvxy(ix_new, iy_new, iz) = temp_qvy(i_old, iy_new, iz)
        endif
      enddo
    enddo
  enddo

  ! Z�������ĵ��ֵ
  do iy_new = 1, ny*nsub + 1
    do ix_new = 1, nx*nsub
      do iz_new = 1, nz*nsub
        z_new = real(iz_new - 0.5) / real(nsub)
        k_old = floor(z_new + 0.5)
        if (k_old < 1) k_old = 1
        if (k_old > nz) k_old = nz
        delta = z_new - (k_old - 0.5)
        if (k_old < nz) then
          v_new(ix_new, iy_new, iz_new) = temp_xy(ix_new, iy_new, k_old)*(1.0 - delta) + temp_xy(ix_new, iy_new, k_old+1)*delta
          qv_new(ix_new, iy_new, iz_new) = temp_qvxy(ix_new, iy_new, k_old)*(1.0 - delta) + temp_qvxy(ix_new, iy_new, k_old+1)*delta
        else
          v_new(ix_new, iy_new, iz_new) = temp_xy(ix_new, iy_new, k_old)
          qv_new(ix_new, iy_new, iz_new) = temp_qvxy(ix_new, iy_new, k_old)
        endif
      enddo
    enddo
  enddo

  deallocate(temp_y, temp_xy, temp_qvy, temp_qvxy)
  return
end subroutine interpolate_vqv !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

    

!############################################################
!                                                           #
! Interpolate W                                             #
!                                                           #
!############################################################ 
subroutine interpolate_w(w, w_new, nx, ny, nz, nsub)
  implicit none
  integer, intent(in) :: nx, ny, nz
  integer, intent(in) :: nsub
  real, intent(in) :: w(nx, ny, nz+1)
  real, intent(out) :: w_new(nx*nsub, ny*nsub, nz*nsub + 1)
  real, allocatable :: temp_z(:,:,:), temp_xz(:,:,:)
  integer :: ix, iy, iz, ix_new, iy_new, iz_new, i_old, j_old, k_old
  real :: x_new, y_new, z_new, delta

  allocate(temp_z(nx, ny, nz*nsub + 1))
  allocate(temp_xz(nx*nsub, ny, nz*nsub + 1))

  ! Z����߽����ֵ
  do iy = 1, ny
    do ix = 1, nx
      do iz_new = 1, nz*nsub + 1
        z_new = real(iz_new - 1) / real(nsub)
        k_old = floor(z_new) + 1
        if (k_old < 1) k_old = 1
        if (k_old > nz) k_old = nz
        delta = z_new - (k_old - 1)
        if (k_old <= nz) then
          temp_z(ix, iy, iz_new) = w(ix, iy, k_old)*(1.0 - delta) + w(ix, iy, k_old+1)*delta
        else
          temp_z(ix, iy, iz_new) = w(ix, iy, k_old)
        endif
      enddo
    enddo
  enddo

  ! X�������ĵ��ֵ
  do iz_new = 1, nz*nsub + 1
    do iy = 1, ny
      do ix_new = 1, nx*nsub
        x_new = real(ix_new - 0.5) / real(nsub)
        i_old = floor(x_new + 0.5)
        if (i_old < 1) i_old = 1
        if (i_old > nx) i_old = nx
        delta = x_new - (i_old - 0.5)
        if (i_old < nx) then
          temp_xz(ix_new, iy, iz_new) = temp_z(i_old, iy, iz_new)*(1.0 - delta) + temp_z(i_old+1, iy, iz_new)*delta
        else
          temp_xz(ix_new, iy, iz_new) = temp_z(i_old, iy, iz_new)
        endif
      enddo
    enddo
  enddo

  ! Y�������ĵ��ֵ
  do iz_new = 1, nz*nsub + 1
    do ix_new = 1, nx*nsub
      do iy_new = 1, ny*nsub
        y_new = real(iy_new - 0.5) / real(nsub)
        j_old = floor(y_new + 0.5)
        if (j_old < 1) j_old = 1
        if (j_old > ny) j_old = ny
        delta = y_new - (j_old - 0.5)
        if (j_old < ny) then
          w_new(ix_new, iy_new, iz_new) = temp_xz(ix_new, j_old, iz_new)*(1.0 - delta) + temp_xz(ix_new, j_old+1, iz_new)*delta
        else
          w_new(ix_new, iy_new, iz_new) = temp_xz(ix_new, j_old, iz_new)
        endif
      enddo
    enddo
  enddo

  deallocate(temp_z, temp_xz)
  return
end subroutine interpolate_w !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@



!############################################################
!                                                           #
! Interpolate W and Qv                                      #
!                                                           #
!############################################################ 
subroutine interpolate_wqv(w, qv, w_new, qv_new, nx, ny, nz, nsub)
  implicit none
  integer, intent(in) :: nx, ny, nz
  integer, intent(in) :: nsub
  real, intent(in) :: w(nx, ny, nz+1)
  real, intent(in) :: qv(nx, ny, nz+1)
  real, intent(out) :: w_new(nx*nsub, ny*nsub, nz*nsub + 1)
  real, intent(out) :: qv_new(nx*nsub, ny*nsub, nz*nsub + 1)
  real, allocatable :: temp_z(:,:,:), temp_xz(:,:,:), temp_qvz(:,:,:), temp_qvxz(:,:,:)
  integer :: ix, iy, iz, ix_new, iy_new, iz_new, i_old, j_old, k_old
  real :: x_new, y_new, z_new, delta

  allocate(temp_z(nx, ny, nz*nsub + 1))
  allocate(temp_xz(nx*nsub, ny, nz*nsub + 1))
  allocate(temp_qvz(nx, ny, nz*nsub + 1))
  allocate(temp_qvxz(nx*nsub, ny, nz*nsub + 1))

  ! Z����߽����ֵ
  do iy = 1, ny
    do ix = 1, nx
      do iz_new = 1, nz*nsub + 1
        z_new = real(iz_new - 1) / real(nsub)
        k_old = floor(z_new) + 1
        if (k_old < 1) k_old = 1
        if (k_old > nz) k_old = nz
        delta = z_new - (k_old - 1)
        if (k_old <= nz) then
          temp_z(ix, iy, iz_new) = w(ix, iy, k_old)*(1.0 - delta) + w(ix, iy, k_old+1)*delta
          temp_qvz(ix, iy, iz_new) = qv(ix, iy, k_old)*(1.0 - delta) + qv(ix, iy, k_old+1)*delta
        else
          temp_z(ix, iy, iz_new) = w(ix, iy, k_old)
          temp_qvz(ix, iy, iz_new) = qv(ix, iy, k_old)
        endif
      enddo
    enddo
  enddo

  ! X�������ĵ��ֵ
  do iz_new = 1, nz*nsub + 1
    do iy = 1, ny
      do ix_new = 1, nx*nsub
        x_new = real(ix_new - 0.5) / real(nsub)
        i_old = floor(x_new + 0.5)
        if (i_old < 1) i_old = 1
        if (i_old > nx) i_old = nx
        delta = x_new - (i_old - 0.5)
        if (i_old < nx) then
          temp_xz(ix_new, iy, iz_new) = temp_z(i_old, iy, iz_new)*(1.0 - delta) + temp_z(i_old+1, iy, iz_new)*delta
          temp_qvxz(ix_new, iy, iz_new) = temp_qvz(i_old, iy, iz_new)*(1.0 - delta) + temp_qvz(i_old+1, iy, iz_new)*delta
        else
          temp_xz(ix_new, iy, iz_new) = temp_z(i_old, iy, iz_new)
          temp_qvxz(ix_new, iy, iz_new) = temp_qvz(i_old, iy, iz_new)
        endif
      enddo
    enddo
  enddo

  ! Y�������ĵ��ֵ
  do iz_new = 1, nz*nsub + 1
    do ix_new = 1, nx*nsub
      do iy_new = 1, ny*nsub
        y_new = real(iy_new - 0.5) / real(nsub)
        j_old = floor(y_new + 0.5)
        if (j_old < 1) j_old = 1
        if (j_old > ny) j_old = ny
        delta = y_new - (j_old - 0.5)
        if (j_old < ny) then
          w_new(ix_new, iy_new, iz_new) = temp_xz(ix_new, j_old, iz_new)*(1.0 - delta) + temp_xz(ix_new, j_old+1, iz_new)*delta
          qv_new(ix_new, iy_new, iz_new) = temp_qvxz(ix_new, j_old, iz_new)*(1.0 - delta) + temp_qvxz(ix_new, j_old+1, iz_new)*delta
        else
          w_new(ix_new, iy_new, iz_new) = temp_xz(ix_new, j_old, iz_new)
          qv_new(ix_new, iy_new, iz_new) = temp_qvxz(ix_new, j_old, iz_new)
        endif
      enddo
    enddo
  enddo

  deallocate(temp_z, temp_xz, temp_qvz, temp_qvxz)
  return
end subroutine interpolate_wqv !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@



!############################################################
!                                                           #
! Interpolate Qvs                                           #
!                                                           #
!############################################################ 
subroutine interpolate_qvs(qv, qv_new, nx, ny, nz, nsub)
  implicit none
  integer, intent(in) :: nx, ny, nz
  integer, intent(in) :: nsub
  real, intent(in)    :: qv(nx, ny, nz)       ! ԭqv���������������ģ�
  real, intent(out)   :: qv_new(nsub*nx, nsub*ny, nsub*nz)  ! ��ֵ���qv���������������ģ�
  integer :: i_parent, j_parent, k_parent     ! ����������
  integer :: ix_sub, iy_sub, iz_sub           ! �������ڸ������ڵľֲ�������0~8��
  integer :: i0, j0, k0, i1, j1, k1           ! ������ǵ�����
  real    :: dx, dy, dz                       ! �����������ڸ������ڵ����λ��
  real    :: wx0, wx1, wy0, wy1, wz0, wz1     ! ��ֵȨ��

  ! ����ÿ��������Ԫ
  do k_parent = 1, nz
    do j_parent = 1, ny
      do i_parent = 1, nx
        ! ����������Ԫ�ڵ�9x9x9��������
        do iz_sub = 0, nsub-1
          do iy_sub = 0, nsub-1
            do ix_sub = 0, nsub-1
              ! --- ����1: ���㸸����ǵ����� ---
              ! x����ǵ�
              i0 = i_parent
              if (i_parent < nx) then
                i1 = i_parent + 1
                dx = (real(ix_sub) + 0.5) / real(nsub)  ! ���λ�� [0.5/9, 8.5/9]
              else
                i1 = i_parent  ! �߽紦�������һ��������Ԫ���Ҳ൥Ԫ
                dx = 1.0       ! Ȩ��ǿ��Ϊ1.0����ʹ��i0��ֵ��
              endif

              ! y����ǵ�
              j0 = j_parent
              if (j_parent < ny) then
                j1 = j_parent + 1
                dy = (real(iy_sub) + 0.5) / real(nsub)
              else
                j1 = j_parent
                dy = 1.0
              endif

              ! z����ǵ�
              k0 = k_parent
              if (k_parent < nz) then
                k1 = k_parent + 1
                dz = (real(iz_sub) + 0.5) / real(nsub)
              else
                k1 = k_parent
                dz = 1.0
              endif

              ! --- ����2: �����ֵȨ�� ---
              wx0 = 1.0 - dx
              wx1 = dx
              wy0 = 1.0 - dy
              wy1 = dy
              wz0 = 1.0 - dz
              wz1 = dz

              ! --- ����3: ��ά���Բ�ֵ ---
              qv_new(nsub*(i_parent-1) + ix_sub + 1, &
                     nsub*(j_parent-1) + iy_sub + 1, &
                     nsub*(k_parent-1) + iz_sub + 1) = &
                wx0 * wy0 * wz0 * qv(i0, j0, k0) + &  ! (i0, j0, k0)
                wx0 * wy0 * wz1 * qv(i0, j0, k1) + &  ! (i0, j0, k1)
                wx0 * wy1 * wz0 * qv(i0, j1, k0) + &  ! (i0, j1, k0)
                wx0 * wy1 * wz1 * qv(i0, j1, k1) + &  ! (i0, j1, k1)
                wx1 * wy0 * wz0 * qv(i1, j0, k0) + &  ! (i1, j0, k0)
                wx1 * wy0 * wz1 * qv(i1, j0, k1) + &  ! (i1, j0, k1)
                wx1 * wy1 * wz0 * qv(i1, j1, k0) + &  ! (i1, j1, k0)
                wx1 * wy1 * wz1 * qv(i1, j1, k1)      ! (i1, j1, k1)
            enddo
          enddo
        enddo
      enddo
    enddo
  enddo
end subroutine interpolate_qvs !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@




!############################################################
!                                                           #
! Operate Subroutine                                        #
!                                                           #
!############################################################ 
SUBROUTINE avgx(a,onvf,nx,ny,nz,ibgn,iend,jbgn,jend,kbgn,kend,aavg)
    implicit none

    integer, intent(in) :: nx, ny, nz                    ! Number of grid points in 3 directions
    real, intent(in) :: a(nx,ny,nz)                      ! Input array
    integer, intent(in) :: onvf                          ! Integer grid point indicator
    integer, intent(in) :: ibgn,iend,jbgn,jend,kbgn,kend ! Integer indicator of multiplication
    real, intent(out) :: aavg(nx,ny,nz)                  ! Averaged array
    
    ! Local variables:
    integer :: i,j,k,iright,ileft
    
    if ( onvf == 1) then
        iright = 0
        ileft = -1
      else
        iright = 1
        ileft = 0
      end if

      do k=kbgn,kend
        do j=jbgn,jend
          do i=ibgn,iend
            aavg(i,j,k)=(a(i+iright,j,k) + a(i+ileft ,j,k))*0.5
          end do
        end do
      end do

  return
END SUBROUTINE avgx


SUBROUTINE avgy(a,onvf,nx,ny,nz,ibgn,iend,jbgn,jend,kbgn,kend,aavg)
    implicit none

    integer, intent(in) :: nx, ny, nz                    ! Number of grid points in 3 directions
    real, intent(in) :: a(nx,ny,nz)                      ! Input array
    integer, intent(in) :: onvf                          ! Integer grid point indicator
    integer, intent(in) :: ibgn,iend,jbgn,jend,kbgn,kend ! Integer indicator of multiplication
    real, intent(out) :: aavg(nx,ny,nz)                  ! Averaged array
    
    ! Local variables:
    integer :: i,j,k,jright,jleft
    
    if ( onvf == 1) then
        jright = 0
        jleft = -1
      else
        jright = 1
        jleft = 0
      end if

      do k=kbgn,kend
        do j=jbgn,jend
          do i=ibgn,iend
            aavg(i,j,k)=(a(i,j+jright,k) + a(i,j+jleft,k))*0.5
          end do
        end do
      end do

  return
END SUBROUTINE avgy



SUBROUTINE avgz(a,onvf,nx,ny,nz,ibgn,iend,jbgn,jend,kbgn,kend,aavg)
    implicit none

    integer, intent(in) :: nx, ny, nz                    ! Number of grid points in 3 directions
    real, intent(in) :: a(nx,ny,nz)                      ! Input array
    integer, intent(in) :: onvf                          ! Integer grid point indicator
    integer, intent(in) :: ibgn,iend,jbgn,jend,kbgn,kend ! Integer indicator of multiplication
    real, intent(out) :: aavg(nx,ny,nz)                  ! Averaged array
    
    ! Local variables:
    integer :: i,j,k,kup,kdown
    
    if ( onvf == 1) then
        kup = 0
        kdown = -1
      else
        kup = 1
        kdown = 0
      end if

      do k=kbgn,kend
        do j=jbgn,jend
          do i=ibgn,iend
            aavg(i,j,k)=(a(i,j,k+kup) + a(i,j,k+kdown))*0.5
          end do
        end do
      end do

  return
END SUBROUTINE avgz
    
    
SUBROUTINE avgz_1d(a,onvf,nz,kbgn,kend,aavg)
    implicit none

    integer, intent(in) :: nz               ! Number of grid points in 3 directions
    real, intent(in) :: a(nz)               ! Input array
    integer, intent(in) :: onvf             ! Integer grid point indicator
    integer, intent(in) :: kbgn,kend        ! Integer indicator of multiplication
    real, intent(out) :: aavg(nz)           ! Averaged array
    
    ! Local variables:
    integer :: k,kup,kdown
    
    if ( onvf == 1) then
        kup = 0
        kdown = -1
      else
        kup = 1
        kdown = 0
      end if

      do k=kbgn,kend
        aavg(k)=(a(k+kup) + a(k+kdown))*0.5
      end do

  return
END SUBROUTINE avgz_1d

